/*
 * =====================================================================================
 *
 *       Filename:  isocontour.cpp
 *
 *    Description:  Generate stroke image from matrix using isocontours and
 *                  morphological thinning
 *
 *        Version:  1.0
 *        Created:  18/07/12 14:18:53
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:  Moos Hueting (mh), mooshueting@gmail.com
 *   Organization:  
 *
 * =====================================================================================
 */
#include <iostream>
#include <cv.h>
#include <highgui.h>
#include <time.h>
#include <math.h>
#include "thin.h"

#define FACTOR 8
#define N_DILATE 1

using namespace std;
using namespace cv;

Mat Skeletonization2(Mat *src) {
    Mat skel(src->size(), CV_8UC1, Scalar(0));
    Mat temp(src->size(), CV_8UC1);

    Mat element = getStructuringElement(MORPH_CROSS, Size(3, 3));

    bool done;
    do {
        morphologyEx(*src, temp, cv::MORPH_OPEN, element);
        bitwise_not(temp, temp);
        bitwise_and(*src, temp, temp);
        bitwise_or(skel, temp, skel);
        erode(*src, *src, element);

        done = (cv::norm(*src) == 0);
    } while (!(done));

    return skel;

}

Mat isocontour(Mat A) {
    Size startSize = A.size();

    Mat mask(startSize.height, startSize.width, CV_64FC1);

    for (int i=0; i < startSize.height; i++) {
        for (int j=0; j < startSize.width; j++) {
            double x = -1.0 + ((i*2)/((double)startSize.width - 1));
            double y = -1.0 + ((j*2)/((double)startSize.height - 1));
            if (pow(x, 4) > pow(y, 4)) {
                if (1-pow(x,4) > 0)
                    mask.at<double>(i,j) = 1 - pow(x,4);
                else
                    mask.at<double>(i,j) = 0;
            } else {
                if (1-pow(y,4) > 0)
                    mask.at<double>(i,j) = 1 - pow(y,4);
                else
                    mask.at<double>(i,j) = 0;
            }
        }
    }

        Mat B(startSize.height, startSize.width, CV_64FC1);
//        GaussianBlur(A, B, Size(7, 7), 1.5);
        A.copyTo(B);

# if 1
        multiply(B, mask, B);
# endif 

        resize(B, B, Size(0,0), FACTOR, FACTOR, INTER_LANCZOS4);

        Size s = B.size();

        for (int i=0; i < s.height; i++) {
            for (int j=0; j < s.width; j++) {
                if (B.at<double>(i, j) > 0.4 && B.at<double>(i,j) < 0.9)
                    B.at<double>(i, j) = 1.0;
                else
                    B.at<double>(i, j) = 0.0;
            }
        }

        B.convertTo(B, CV_8UC1);
        B *= 255;

        B = thin(B);

        //Mat strel = getStructuringElement(MORPH_ELLIPSE, Size(3,3));

        //dilate(B, B, strel, Point(-1,-1), N_DILATE);
        //B = thin(B);
        //dilate(B, B, strel, Point(-1,-1), N_DILATE);

        B.convertTo(B, CV_64FC1);
        B = B / 255.0;

        return B;
}
